{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "9fa20760",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "num=1001\n",
    "length=10\n",
    "x = np.random.random(size=num)*length  #分子位置随机化0-10之间\n",
    "y = np.random.random(size=num)*length\n",
    "v=np.zeros(num)+2.0  #初速度为1，方向0-360平均分布\n",
    "o=np.arange(2*np.pi/num,2*(num+1)*np.pi/num,2*np.pi/num)\n",
    "vx=v*np.cos(o)\n",
    "vy=v*np.sin(o)\n",
    "dt=0.1\n",
    "radius=0.2 #分子见距离小于radius认为发生碰撞\n",
    "\n",
    "def molecularcollision(x1,y1,x2,y2,vx1,vy1,vx2,vy2):  #计算两分子碰撞后的速度\n",
    "    #转换坐标系到质心连线为x' 的速度v1temp,v2temp\n",
    "    angle=np.arctan(abs(y1-y2)/abs(x1-x2))\n",
    "    vx1Temp=vx1*np.cos(angle)+vy1*np.sin(angle)\n",
    "    vy1Temp=vy1*np.cos(angle)+vx1*np.cos(angle+np.pi/2)\n",
    "    vx2Temp=vx2*np.cos(angle)+vy2*np.sin(angle)\n",
    "    vy2Temp=vy2*np.cos(angle)+vx2*np.cos(angle+np.pi/2)\n",
    "    #互换并转换为原坐标系\n",
    "    \n",
    "    vx1=vx2Temp*np.cos(angle)+vy1Temp*np.cos(angle+np.pi/2)\n",
    "    vy1=vx2Temp*np.sin(angle)+vy1Temp*np.sin(angle+np.pi/2)\n",
    "    vx2=vx1Temp*np.cos(angle)+vy2Temp*np.cos(angle+np.pi/2)\n",
    "    vy2=vx1Temp*np.sin(angle)+vy2Temp*np.sin(angle+np.pi/2)  \n",
    "    \n",
    "    return vx1,vy1,vx2,vy2\n",
    "    \n",
    "def collisionDetect(x,y,vx,vy,dt,radius,length): #检测是否发生碰撞，碰撞有边界碰撞，分子间碰撞\n",
    "    #边界碰撞，现在位置x,y计算dt后的x,y,当x，y超出边界，超出的反向\n",
    "    #边框边界 0<x<10,0<y<10,x超x反向，y超y反向\n",
    "    xtemp=x+dt*vx\n",
    "    ytemp=y+dt*vy\n",
    "    for i in range(len(x)):\n",
    "        if xtemp[i]>0 and xtemp[i]<length:\n",
    "            vx[i]=vx[i]\n",
    "        else:\n",
    "            vx[i]=-vx[i]\n",
    "        if ytemp[i]>0 and ytemp[i]<length:\n",
    "            vy[i]=vy[i]\n",
    "        else:\n",
    "            vy[i]=-vy[i]\n",
    "            \n",
    "    #增加两条线段，线段1 （3，4） （3，6），线段2 （7，4） （7，6）\n",
    "    #判断分子是否碰到了线段，首先y和ytemp中一个在4-6之间，另外x和xtemp中一个大于3，一个小于3\n",
    "    lCollision=rCollision=0\n",
    "    for i in range(len(x)):\n",
    "        if (y[i]>4*length/10 and y[i]<6*length/10) or (ytemp[i]>4*length/10 and ytemp[i]<6*length/10):\n",
    "            left=3*length/10\n",
    "            if (x[i]<left and xtemp[i]>left) or (x[i]>left and xtemp[i]<left): #向碰撞线段3\n",
    "                vx[i]=-vx[i]\n",
    "                lCollision=lCollision+1\n",
    "            right=7*length/10\n",
    "            if ( x[i]<right and xtemp[i]>right)or (x[i]>right and xtemp[i]<right):#向碰撞线段7\n",
    "                vx[i]=-vx[i]\n",
    "                rCollision=rCollision+1\n",
    "    \n",
    "    #为了简化计算认为 |x1-x2| |y1-y2|同时小于radius就认为发生了碰撞交换\n",
    "    s=0\n",
    "    for j in range(len(x)):\n",
    "        s=s+1\n",
    "        for k in range(s,len(x)):\n",
    "            dx=abs(x[j]-x[k])\n",
    "            dy=abs(y[j]-y[k])\n",
    "            if dx<radius and dy<radius:\n",
    "                vx[j],vy[j],vx[k],vy[k]=molecularcollision(x[j],y[j],x[k],y[k],vx[j],vy[j],vx[k],vy[k])\n",
    "    return vx,vy,lCollision,rCollision\n",
    "\n",
    "def runAndPlot(n,x,y,vx,vy,dt,radius,length):\n",
    "    l=r=0\n",
    "    for i in range(n):\n",
    "        vx,vy,lCollision,rCollision=collisionDetect(x,y,vx,vy,dt,radius,length)\n",
    "        l=lCollision+l\n",
    "        r=rCollision+r\n",
    "        x=x+vx*dt\n",
    "        y=y+vy*dt\n",
    "        ''' \n",
    "        plt.clf()\n",
    "        plt.xlim(0, 10) \n",
    "        plt.ylim(0, 10) \n",
    "        plt.figure(figsize=(6, 6))\n",
    "        \n",
    "        plt.plot([3,3],[4,6])\n",
    "        plt.plot([7,7],[4,6])\n",
    "        \n",
    "        plt.scatter(x, y)\n",
    "        plt.scatter(x[0], y[0],c='r')\n",
    "        plt.savefig('molecule/test'+\"{:0>3d}\".format(i)+'.jpg') '''\n",
    "        if ((i+1) % 200)==0:\n",
    "            print(\"左边线段碰撞次数\",l,\"右边线段碰撞次数\",r)\n",
    "    print(\"左边线段碰撞次数\",l,\"右边线段碰撞次数\",r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "358d8204",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "左边线段碰撞次数 413 右边线段碰撞次数 454\n",
      "左边线段碰撞次数 823 右边线段碰撞次数 888\n",
      "左边线段碰撞次数 1232 右边线段碰撞次数 1316\n",
      "左边线段碰撞次数 1655 右边线段碰撞次数 1751\n",
      "左边线段碰撞次数 2126 右边线段碰撞次数 2160\n",
      "左边线段碰撞次数 2567 右边线段碰撞次数 2578\n",
      "左边线段碰撞次数 2963 右边线段碰撞次数 3002\n",
      "左边线段碰撞次数 3378 右边线段碰撞次数 3432\n",
      "左边线段碰撞次数 3795 右边线段碰撞次数 3826\n",
      "左边线段碰撞次数 4213 右边线段碰撞次数 4284\n",
      "左边线段碰撞次数 4655 右边线段碰撞次数 4701\n",
      "左边线段碰撞次数 5118 右边线段碰撞次数 5139\n",
      "左边线段碰撞次数 5534 右边线段碰撞次数 5600\n",
      "左边线段碰撞次数 5950 右边线段碰撞次数 6016\n",
      "左边线段碰撞次数 6372 右边线段碰撞次数 6409\n",
      "左边线段碰撞次数 6800 右边线段碰撞次数 6828\n",
      "左边线段碰撞次数 7223 右边线段碰撞次数 7257\n",
      "左边线段碰撞次数 7661 右边线段碰撞次数 7652\n",
      "左边线段碰撞次数 8087 右边线段碰撞次数 8076\n",
      "左边线段碰撞次数 8509 右边线段碰撞次数 8519\n",
      "左边线段碰撞次数 8944 右边线段碰撞次数 8968\n",
      "左边线段碰撞次数 9336 右边线段碰撞次数 9404\n",
      "左边线段碰撞次数 9771 右边线段碰撞次数 9796\n",
      "左边线段碰撞次数 10202 右边线段碰撞次数 10212\n",
      "左边线段碰撞次数 10621 右边线段碰撞次数 10637\n",
      "左边线段碰撞次数 11061 右边线段碰撞次数 11051\n",
      "左边线段碰撞次数 11446 右边线段碰撞次数 11504\n",
      "左边线段碰撞次数 11855 右边线段碰撞次数 11935\n",
      "左边线段碰撞次数 12294 右边线段碰撞次数 12350\n",
      "左边线段碰撞次数 12737 右边线段碰撞次数 12785\n",
      "左边线段碰撞次数 13145 右边线段碰撞次数 13231\n",
      "左边线段碰撞次数 13558 右边线段碰撞次数 13632\n",
      "左边线段碰撞次数 14000 右边线段碰撞次数 14034\n",
      "左边线段碰撞次数 14429 右边线段碰撞次数 14456\n",
      "左边线段碰撞次数 14837 右边线段碰撞次数 14884\n",
      "左边线段碰撞次数 15273 右边线段碰撞次数 15306\n",
      "左边线段碰撞次数 15705 右边线段碰撞次数 15737\n",
      "左边线段碰撞次数 16157 右边线段碰撞次数 16155\n",
      "左边线段碰撞次数 16603 右边线段碰撞次数 16580\n",
      "左边线段碰撞次数 17005 右边线段碰撞次数 16989\n",
      "左边线段碰撞次数 17431 右边线段碰撞次数 17411\n",
      "左边线段碰撞次数 17853 右边线段碰撞次数 17844\n",
      "左边线段碰撞次数 18283 右边线段碰撞次数 18264\n",
      "左边线段碰撞次数 18748 右边线段碰撞次数 18695\n",
      "左边线段碰撞次数 19169 右边线段碰撞次数 19116\n",
      "左边线段碰撞次数 19578 右边线段碰撞次数 19542\n",
      "左边线段碰撞次数 19997 右边线段碰撞次数 19965\n",
      "左边线段碰撞次数 20396 右边线段碰撞次数 20392\n",
      "左边线段碰撞次数 20801 右边线段碰撞次数 20808\n",
      "左边线段碰撞次数 21232 右边线段碰撞次数 21220\n",
      "左边线段碰撞次数 21232 右边线段碰撞次数 21220\n"
     ]
    }
   ],
   "source": [
    "runAndPlot(10000,x,y,vx,vy,dt,radius,length)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "84031d17",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
